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ABSTRACT 


In double differencing a regression system obtained from concurrent Global 
Positioning System (GPS) observation sequences, one either under-samples the 
system to avoid introducing colored measurement statistics, or one fully samples 
the system incurring the resulting non-diagonal covariance matrix for the 
differenced measurement errors. A suboptimal estimation result will be obtained 
in the under- sampling case and will also be obtained in the fully sampled case 
unless the color noise statistics are taken into account. The latter approach 
requires a least squares weighting matrix derived from inversion of a 
non-diagonal covariance matrix for the differenced measurement errors instead of 
inversion of the customary diagonal one associated with white noise processes. 
This publication presents the so-called fully redundant double differencing 
algorithm for generating a weighted double differenced regression system that 
yields equivalent estimation results, but features for certain cases a diagonal 
weighting matrix even though the differenced measurement error statistics are 
highly colored. 
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INTRODUCTION 


The concurrent operation of two or more Global Positioning System (GPS) ground 
terminals observing multiple GPS satellites generates multiple data streams that 
yield high accuracy relative positioning. Concurrent operations produce 
simultaneous observational epoch sequences among the multiple data streams, 
which allows one to eliminate errors in the state vector estimates that arise 
from common mode-error sources in the data streams. In particular, the effects 
of unknown and imperfectly running clocks on board the satellites and in the 
ground terminals may be avoided. Although a parameterized representation of the 
clock errors in the regression equations may be estimated in parallel with the 
state vector in an expanded state space mode, this approach leads to very large 
arrays unless a partitioning technique is used. Alternative approaches 
eliminate the clock terms before obtaining the estimate of the state vector. In 
the latter approach, the regression equations are double differenced at 
simultaneous epochs to cancel the clock parameters before the state vector is 
estimated. Alternatively, the regression equations are not differenced; rather, 
a series of orthogonal transformations are applied to the system. These 
transformations effectively eliminate these common mode errors while yielding a 
minimum variance estimate of the state vector. One technique following the 
latter approach involves the operation on the regression system in either a 
batch or sequential mode by a series of Householder matrices, tailored to this 
problem, that transforms the information matrix into the form of an upper 
triangular array. The lower portion of this upper array is explicitly free of 
clock information. Thus, these Householder transformations effectively 
eliminate clock errors before the actual inversion of the state vector 
triangularized information matrix takes place [1,2]. The advantage of the 
Householder matrices or similar approaches is that additional modeling and a 
priori information about these errors can be readily introduced in the 
estimation process without explicitly solving for the clock errors. For 
example, the stochastic properties of the clock variability may be modeled by 
specification of the parameters from a first order Markov process. 

Alternatively, the clocks in GPS terminals collocated at VLBI fiducial sites can 
be slaved to the resident hydrogen maser. Their near-linear time behavior 
provides important additional information, particularly in the control of errors 
in the estimated GPS ephemerides, which would be lost in a double differencing 
mode. 
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The double differencing approach appears to offer simplicity by eliminating the 
clock parameters straightaway. Double differencing is attractive because the 
differenced residuals lend themselves to an easy and traceable physical 
interpretation, which is not usually the case with residuals obtained from 
Householder or similar orthogonal transformations. Double differencing is also 
particularly useful in a noisy environment for editing functions, such as 
detecting and removing carrier cycle dropouts and other error sources from the 
data streams. A potential disadvantage of double differencing is that it 
usually results in colored noise in the estimation problem. Also, techniques 
that involve an orthogonalization of the information matrix such as, for 
example, the Householder approach, tend to be more computationally efficient. 
Nevertheless, there are situations where double differencing may be a practical 
necessity or even preferable in a data management sense. 

This publication addresses two full-rank approaches to double differencing: the 
so-called minimally redundant approach and the fully redundant approach. A 
full-rank approach is one in which the double differencing in a matrix operator 
sense has maximal rank and hence spans the linear vector space defined by the 
regression system. Double differencing introduces colored noise statistics into 
the differenced measurement error covariance. To obtain estimation results that 
are equivalent to those obtained by the expanded state or Householder methods, 
these approaches require a least-squares weighting matrix that is derived either 
from an inversion of the colored-measurement error-covariance matrix, which is 
nondiagonal to a varying degree depending on the differencing scheme used and 
the number of stations and satellites involved. Alternatively, one can use a 
whitening transformation on the differenced regression system before obtaining 
the inverse. 

The computational time for numerical inversion of the nondiagonal covariance 
matrix tends to grow as the cube of its size. However, if the undifferenced 
measurement statistics are white, then the double-differenced covariance matrix 
may be only mildly striped if either the number of terminals or satellites is 
moderate. In this case, more efficient matrix inversion techniques are 
available [3], which tend to grow somewhat more slowly than the cube of the 
matrix size. Also, analytic inverses are available for certain measurement 
error covariances with simple forms. The fully redundant approach features a 
diagonal weighting matrix for certain forms for the undifferenced 
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measurement-error covariance even though the double-difference measurement-error 
statistics are highly colored. This feature may provide a computational 
advantage in limited situations. On the other hand, the number of differencing 
operations themselves in this approach tends to grow q[uadratically with the 
number of terminals and with the number of GPS satellites. Comparative studies 
indicate that this method should be considered competitive in computational time 
only in those cases where the number of terminals or the number of observed 
satellites is modest, i.e., roughly four or less, and only in cases where 
analytic forms for the diagonal weighting matrix are available. If analytical 
forms are not available and one wishes to use double differencing, then a 
full-rank but minimally redundant doubled-differencing approach appears to be 
the most straightforward [3]. 

After a general description of the regression problem given in the next section, 
this publication discusses double differencing operations and the equivalence 
under certain conditions of the resulting least-squares estimates to those 
obtained by other techniques. The publication then discusses minimally and 
fully redundant double difference matrix operators and their least-squares 
weighting matrices. In the fully redundant approach, its corresponding 
white-noise-equivalent diagonal weighting matrix is developed. It is shown that 
analytic forms of this .diagonal matrix can be generated in some cases. In 
theory a composite diagonal and antisymmetric weighting matrix can be generated 
in all cases, but the practicability of generating analytic forms appears 
limited to certain restrictive assumptions about the topological properties of 
the undifferenced measurement error covariance. Finally, the paper treats the 
problem of missing data streams or outages in double differencing. 
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II. PROBLEM STATEMENT 


We consider the case where m ground terminals are concurrently tracking r GPS 
satellites from which mr distinct data streams are generated in parallel with 
simultaneous observational epoch sequences in each stream. The measurement 
types assumed here are either the undifferenced carrier phase or the group delay 
based ranging, which are characteristic of single GPS receiver operations. The 
carrier phase measurements are assumed to be phase connected although this 
assumption is not fundamental to what follows. In addition to depending on the 
state vector, the measurements are corrupted by the offsets of the clocks from 
universal time. A linear regression system of the form 

y + e = [B:H] 


+ e 


( 1 ) 


is assumed where the vectors y and X represent small differences from nominal 
values and where: 

y is the mrN by 1 observation vector 

N is the number of observational epochs in each data stream 
^ is the grand information matrix for the regression system 
X is the expanded state vector of dimension p + (m + r)N by 1 
e is the measurement error vector 

B is the mrN by (m + r)N clock vector information matrix 
H is the mrN by p state vector information matrix 
b is the (m + r)N by 1 clock error parameter vector 
X is the p by 1 state vector 

The state vector consists of the parameter set except clocks that are required 
to properly define the dynamical and observational systems. Thus, in addition to 
state vector information for GPS terminals and GPS satellites, x would also 
include carrier cycle ambiguities, propagation media parameters, etc. 
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The clock characteristics aboard each satellite and in each terminal are not 

modeled in this problem. Rather, each clock error at any epoch is assumed 

(rather pessimistically) to be uncorrelated with any other. Thus, the vector b 

represents the (m + r)N separate clock errors. The rank ofBis(m+r-l)N 

because the sum of the clock parameters at each epoch is unobservable for the 

linear problem. The information matrix H is assumed to be of rank p. The error 

vector e is assumed to have standard mean zero white Gaussian noise properties 

and a covariance matrix of A . 

e 


It is assumed that y and tx/ have been appropriately formatted, with the 
time-ordered subblock of observation residuals and regression coefficients 
generated from observations of the same satellite by the same terminal grouped 
in contiguous rows. Additionally, the m subblocks arising from observations of 
the same satellite by the m different terminals are placed contiguously to 
create a single block containing the information from the entire set of 
observations of the same satellite. The r blocks associated with observations 
of the r different satellites, each with the identical substructure, are grouped 
contiguously to fill out the complete observation vector and information matrix. 
For example, with r = 4 and m = 3, H would appear as 
m = 3 


T l T** T'* rp» rp*« rp« rp« rp 

I ^ M ^ ^ ^ A ^ ^ ^ JL 9 ^ 9 ^ 9 9 M 9 w 9 M 9 9 JL 9 

[^ii*2i^3ii®il®2l®3ii^ii^2iSii°ii°2i°3 


( 2 ) 


r = 4 

where the subscript denotes the terminal and the letters A, B, C and D denote 
the information submatrices from N observations each from four different 
satellites, denoting subscripts as terminal designations and superscripts as 
satellite designations, then, following the form in Eg. (2), H may be written as 


m subblocks 



( 3 ) 


r blocks 
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The units and scale of the clock parameters have been chosen the same as those 
of the observation vector y. It follows, provided that the regression equations 
have been formatted according to Eg. (3), that B, the clock parameter 
information matrix, consists of ones and zeros in the pattern shown below 



The identity matrices in Eqs. (5) and (6) are N by N. 

In a batch mode version of the Householder approach for eliminating the clock 
parameters [1], one operates on Eg. (1) with a series of orthogonal Householder 
matrices [2] that transforms it into the form 
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where R is a square (m+r+p)N matrix consisting of an upper triangular array with 
zeros below the diagonal. In particular, the minimum variance state vector 
estimate x is given by 


X = y 

X "^x 


( 8 ) 


and its covariance matrix A is given by 


= R^ (r'^)'^ 

X X X 


(9) 


It is noted that R is the transformed information matrix in which the clock 

X 

information is explicitly absent; thus, there is no need to jointly solve for 
clock parameters unless one desires their estimates or has other a priori 
information about them. 


A weighted least squares method that is equivalent to but less efficient than 
the Householder method is the so-called augmented state space approach. Because 
it is relatively better known we will use this formulation rather than the 
Householder formulation to establish the equivalency of the double differencing 
approach. Here, one adjoins the clock parameter vector b, to the state vector x, 
and obtains a joint minimum variance estimate of this augmented state vector 
using a standard least squares batch mode approach. The grand covariance matrix 
of the augmented state vector estimate is obtained from the pseudoinverse of 



e 


h^a’ 

e 


I 
I 

I 

a’’ H ! B^ A"’ 
e 1 e 


( 10 ) 


The pseudoinverse of Eg. (10) can be written in partition form as 



Aa 

X 


F 


G 


( 11 ) 
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where is the covariance of the state vector estimate x, and where F and G are 

the cross-covariance and covariance associated with the clock parameter 

T 

sequences. (Pseudoinverses appear because the rank of B B is m+r-1 instead of 

m+r. By fixing the clock at the first ground terminal, for example, one then 

T 

eliminates the first column from B, thereby rendering B B into a full rank 
matrix. In this case, the pseudoinverse operations in the subsequent discussions 
may be replaced by a standard inverse. The ensuing results are not materially 
different from those obtained with the pseudoinverse approach.) We assume that 
H is of full parameter rank p and that the column vectors of H are linearly 
independent of those in B. Using the properties of pseudomatrices, it can be 
shown for matrices of the form in Eq. (11) that 

A- [h’^hh] - d2) 

where the matrix II is defined by 


n = - A"^ bib'^a'^b^ b'^a"^ 

e e ^ e " e 


(13) 


Hence, if the inverse of [H IlH] exists, it follows that covariance of x is given 
by 


= h'^iih 

X 


Also, it follows that x is given by 


(14) 


^ —IT 

X = A- H^ny 


(15) 


It can be shown [2] that Eqs.(8) and (15) are equivalent and that Eqs.(9) and 
(14) are also equivalent. Even in this augmented state space approach, by using 
this matrix partitioning technique one need not explicitly solve this system for 
the clock parameters. 


Because the matrix II and its factorization serves as a bridge between the 
augmented state space and the double differencing approaches, a few observations 
about some of its properties should be noted. From the least squares projection 
theorem in linear vector spaces [4] (or from multiplying Eq. (13) by B), it 
follows that n is orthogonal to B. Also, II, which has a rank of (m-l)(r-l)N, has 
the property that 
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( 16 ) 


n =n A^n 

e 

The matrix II , given by 

n = aJh (17) 

is an idempotent matrix and its own pseudoinverse with a trace equal to 
(m-l)(r-l)N. 

Finally, it should be noted that II itself, which is a relatively simple matrix 
to generate, serves as a matrix operator which, when applied to Eg. (1), 
transforms it into a new regression system that is explicitly free of cloclc 
terms. By premultiplying Eg. (1) by II , one obtains 

ny=IlHx+IIe (18) 

It is easily shown that the resulting least sguares estimate for x and its 
covariance matrix obtained from Eg. (18) are formally identical to those in Egs. 
(15) and (14). This suggests yet another method in lieu of double differencing 
for eliminating clock parameters. 


III. DOUBLE DIFFERENCING OPERATIONS AND CONDITIONS FOR EQUIVALENCE 

In the double differencing approach one operates on the system in Eg. (1) in 
such a way to explicitly eliminate the clock parameters from the resulting 
regression system before carrying out the filtering operation. 

A double differencing operation first differences the regression eguations 
associated with the simultaneous observations by two terminals of the same 
satellite, and then it differences two such first differences obtained from the 
same two terminals simultaneously observing two different satellites. 
Alternatively, these operations may be reversed, but the result is the same. 
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While double differencing operations can be variously formulated, all 
formulations have the property of totally eliminating the clock matrix B from 
the resulting double differenced regression equations. In matrix notation, let 
be a double difference matrix operator. Then, it follows that 

SB = 0 (19) 

If one carries out the double differencing of the regression system in Eq. 
(1) following theS-matrix formulation, one obtains 

2)y =SHx +Se (20) 

a new system of regression equations to be solved in a least squares sense for 
the state vector x. Letting e be defined by 


e=Se 


( 21 ) 


we have for the covariance of e 


A - Sas'’’ (22) 

^ e 

e 

a nondiagonal, colored matrix whose size depends on one*s choice for The 

rank of is at most mrN, which may be smaller than its size depending on the 
e 

specific differencing algorithm used. Consequently, one should use the 
pseudoinverse in the formulation for the least squares estimate of x. 

Then, A - is given by 


and X is given by 


X L e J 




(23) 


(24) 


If the state vector estimate and its covariance in Eqs. (24) and (23) are to be 
equivalent to those in Eqs. (15) and (14), it follows that 
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( 25 ) 


To establish the conditions for this equivalence, it is sufficient to find mrN 
linearly independent nonzero vectors ^ such that 

Mt = 0 (26) 

The matrix B contains (m + r 1)N linearly independent column vectors and from 
the orthogonality condition in Eq. (19) and the properties of pseudoinverse 
matrices it follows that 

MB = 0 (27) 

Next, we premultiply Eq. (25) ^ ^ (assumed to be full rank) and then by 3) to 

obtain 

SJA^M 5 0 (28) 

If 3) fully spans the vector space orthogonal to the space spanned by the columns 
of B, it will contain (m - l)(r - 1)N linearly independent column vectors. 
Because is also orthogonal to B, Eqs. (27) and (28) serve to establish that M 
= 0 and, thus, to establish the equivalence of Eqs. (23) and (14), also Eqs. 

(24) and (15) When’S^ fully spans the space orthogonal to B. It is readily shown 
that any double difference operator that preserves the white noise statistics of 
the differenced measurement errors must have a rank less than (m-l)(r-l)N; 
hence, in general it will yield non-equivalent results. 


An example of a full rank but minimally redundant operator is given by 



(m-l)(r-l) rows 


(29) 


where D, a first difference operator of full rank m-1 but minimally redundant, 
is given by 
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(m-1) rows 


( 30 ) 


1 -1 0 • • • 0 
i 0-1 •• • 0 


1 0 0 • • • -1 


The operator in Eq. (29) is minimally redundant because it generates the least 
number of rows in the double differenced regression matrix while maintaining a 
full rank of (m-l)(r-l). It is called "redundant" because one row and one block 
of the undifferenced system appear respectively in all rows and all blocks of 
the differenced system. The covariance of e (Eq. (22)) using this operator is 
non-diagonal and generally requires a numerical inversion to obtain the least 
squares weighting matrix. However, for modest values of m or r and/or for 
certain simple forms of A^, analytic inverses can be readily derived. (See 
Appendix 1 . ) Remondi [ 3 ] has provided an example of this for m = 2 . 


IV. THE FULLY REDUNDANT DOUBLE DIFFERENCE MATRIX OPERATOR 


Because the double differenced measurement error covariance matrix in Eq. (22) 
may be highly non-diagonal as a result of the differencing operations and 
because it may also be quite large, we present an algorithm which enables one to 
replace it with an equivalent diagonal form. To accomplish this we define a 
specific form for the double difference matrix operator. (See Appendix 2 for a 
more detailed discussion. ) 


Let us define a first difference operator Dj^, which when applied to another 
matrix of k block rows, generates the distinct pair differences that can be 
formed from these block rows. For example, for k = 1, 2, 3, and 4, is given 
by 


D, = 0, Dj =[i ; -1] , D3 = 


I 




I I 0 I -I ' 

0 1 I I -I 






1 t 

1 

-I 1 

0 

0 


— (- 

— 



I 1 

0 1 



-I 

0 

I 1 

n 
0 1 

0 

-I 


--H- 

— 

— 

0 1 

I 1 

-I 

0 

0 1 

1- 

I 1 

0 

-I 

1- 

0 1 
1 

— 1 
0 j 

I 1 

1 

1 -I 

1 

* 





(31) 
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Thus, Dj^ operating on a k by P matrix H creates a new ^ 2 ) ^ matrix, whose 

rows comprise the ^ 2 ) distinct ways that pair differences can be formed from the 
original k rows. 


Next, the analogous double difference operator S) is defined by 



mrN columns 


(32) 


(Henceforth, D will be interpreted as an mth order D-matrix unless otherwise 
specified; thus, the subscript m is dropped.) The rank of is (r-l)(m-l)N; 
hence, fully spans the vector space orthogonal to B and therefore, the 
operation on H by fully retains the information content in H. 

The operation H converts the undifferenced information matrix H, consisting of 
mrN rows, into a fully redundant double differenced information matrix, H, 
consisting of { 2 ) ^ rows. In the example of Eq. (2) with r = 4 and m = 3, 
the operation H yields 18 row blocks of the form 




(A^~ A 2 “ 

(B, - - C, + Cj)! 

(A| A^ “ B| + 
(A 2 -A 3 -B 2 +B/ 


(A, -A2-C, + c/ 

(B, - B2 - D, + Dj)^ 


(A, -A2-D, + D2)^ 1 
(C, - Cj - D, + D2)T i j 

II 

(C2 — C3 — D2 + Dj)^ 1 1 


(33) 


One of the potential disadvantages of this double differencing algorithm for 
large m and r is clearly seen from Eq. (32) or (33); the number of rows 
generated byS> increases quadratically in both m and r. This should be compared 
with the Householder approach which does not expand the regression system. In a 
sense, the ^ operator generates extra or "redundant" 

regression blocks; but, as we shall now show, these lead to a diagonal form that 
is equivalent to 
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V. 


DIAGONAL FORMS FOR THE -MATRIX FORMULATION 


From the equivalence of Eqs. (14) and (23), it follows that II, defined by Eq. 
(13), is also given by 

In the case where A is diagonal and uniform, II can be reduced to a remarkably 

e 

simple form given by 

^T ^ (35 

n=®^w2) 


were W is a 
treating 


(2)(2)^ diagonal weighting matrix. 


as white and replacing 




with W in Eqs. 


This is equivalent to 
(23) and (24). 


This raises the question whether or not other structures also lead to 

diagonal forms. We will show that this is so and we will obtain analytic 
expressions for W for certain cases. 

The Case of Uniform Measure Error Covariance 


For this case, A is given by 
e 

A = 0^1 
e 

where I is an mrN x mrN identity matrix, 
pseudoinverse of is given by 


(36) 


It is shown in Appendix 1 that the 



mr 


T 


(37) 


For this case, II in Eq. (34) reduces to 


TD ra 

which, using the pseudoinverse properties of 21 again, reduces to 

n= — ^ 2)^2 

mra 


(38) 


(39) 
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T T 2 2 2 2 

Thus, because of the rank deficiency of Sf] , both [3>3) ]/m r a and l/mro 

^ serve as the pseudoinverse in Eq. (34) and provide the same result for II. 

Factorization of II in Eq. (39) leads to a form for given by 




( 40 ) 


where the weighting matrix W in this case is given by 


W = — ij- I (41) 

mro 

Also, X is given by 

An alternative derivation based on an explicit evaluation of Eq. (13) is 

provided in Appendix 3. It is noted that Eq. (38) provides an interesting 

T T 2 2 2 

duality property. If one thinks of Eq. (38) as the form® [®2> ^ o » then 

one has the double differencing formulation involving the inner products of 
arrays of dimension ^2)(^) other hand, if one thinks of Eq. (38) 

having the form then one has the equivalent matrices that 

must result from the expanded state or from the Householder approaches. Here, 
the array products are only mrN in dimension. 


®hJ*^|^ W^®yj 


(42) 


1 T 

An operational definition of f® 2D is provided by Equation (A-14) in Appendix 

1 T 

2. An inspection of Eq. (A- 14) shows that the operation by — [® ®] on the 
matrices H and y involves the double subtraction of two mean values of the 
elements of these matrices from each element (plus a correction for double 
counting). These mean values are obtained from averaging over the set of 
terminals and over the set of satellites. Bierman fS] has shown that the 
Householder operation leading to triangularization involves the identical 
averaging processes in the case of uniform measurement error covariance. These 
averaging and subtraction processes have the net effect of cancelling the effect 
of clock terms in y except that it is accomplished by a linear combination of 
the rows of y rather than by explicit double differencing. 
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The Case When A Is Factorable 
e 


When A is factorable into the form 
e 


% 


i 

Q 

0 

1 


- 


b 0 

m 


\ 

\ 

\ 


\ 

\ 


\ 

\ 


\ 

\ 


\ 

0 I 


0 

cr 

m 


J 


(43) 


where « is a scalar and b is an mN by mN diagonal matrix, or equivalently, when 
the standard deviatior 
satellite is given by 


the standard deviation of the error for the jth terminal observing the v th 


V V, (44) 

then it may be shown using techniques described in Appendices 2 and 3 that a 
diagonal form for W is given by 


w - llA^ir^ 


112 2 
0, Oj o, Oj 


112 2 
o, 03 a, 03 


1 12 2 
_ 1 <y_ _ 1 


(J)N, 


113 3 

®1 ®2 ®1 ‘’2 


I 13 3 

I O 


f-1 r-1 rr 
^2 ^2 


r-1 L-1 r r 


i-r 


(45) 


Eq. (45) is also obtained by renormalizing the clock parameters. These 
rescaling factors provide m+r-1 degrees of freedom, which are sufficient to 
transform this problem into an equivalent uniform one. For the factorable case 
it may be shown that II has the form 
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( 46 ) 


II = S!=a“^-a“^ ii''^a“V|| a^II^ 

e e e ■ • e" * 

where 1 is the N by mrN unitary column matrix (see Appendix 1) and 
II Ae^ll Euclidean norm of Eq. (45) is not valid when has a 

general diagonal form. 

The success in deriving a simple analytical expression for the diagonal 
form of W in the cases of uniform and factorable measurement error covariances 
encourages one to seek other cases where this is so. Wu [6] has shown for 
problems involving single differenced observables with one independent clock 
error sequence for each data stream, for example, VLBI-like sequences, that an 
analogous diagonal form for this single differenced version always can be 
established when A ^ is diagonal. In the double differenced case, the matter is 
less straightforward. 

VI. THE EXISTENCE OF DIAGONAL FORMS FOR THE 2J-MATRIX FORMULATION 

To explore the existence question of possible diagonal forms for W one equates 
the forms for IT as given in Eq. (13) and in Eq. (35); from this one derives the 
appropriate form for W. For certain forms for A ^ one can actually carry out 
the matrix algebra operations indicated in Eq. (13) (See Appendix 3). 
Alternatively, one can use the idempotent-related properties of 11 given in Eqs. 
(16) and (17) (see Appendix 4). In the more general case II would be generated 
numerically using Eq. (13); the numerical values of the elements of W could 
still be obtained from equating Eqs. (13) and (35). However, any efficiencies 
gained in using the fully redundant approach would be lost in the numerical 
inversion of Eq. (13). From Eq.(13) II may be written in the form 


r *1 


A 1 B 1 C 1 E 


T T r 


I F G j H 


t' 



1 

1 — 

1 — 

X 

h— 

LU 

1 



where each of these mN by mN sub-matrix blocks of II is constrained by 

IIb'*' = IIB=0 (48) 


17 



2 2 2 2 2 2 

which constitutes (m r N -(m-1) (r-1) ]N independent orthogonality conditions. 

Arbitrary forms for break the symmetries in the off-diagonal sub-blocks of II. 

Although II is symmetric, its individual off-diagonal submatrix blocks are in 

general not symmetric. It would follow that W could not be diagonal. However, 

T 

because of the rank deficiency of [D WD], the form of W is not unique, as was 

T 

already demonstrated in the uniform case. Moreover, inasmuch as H IIH is an 

array of quadratic forms, II can be reformatted without altering the values of 
T 

the elements in H IlH. This is accomplished by expressing each sub-block of II 
as the sum of a symmetric matrix and an antisymmetric one. For example, the 
matrix E in Eq. (47) may be written as 

E = g + E (49) 

where I and E are respectively the symmetric and antisymmetric components of E, 
with 

g = (E + E*^)/2 (50) 

and 

E = (E - E*^)/2 (51) 

With these replacements II may be written as 


n = n + n 


(52) 


where n and II are symmetric. The forms of these matrices are given by 


and 


IT = 


A 

B 

C 

g 


B 

F 

5 

g 


C 

g 

I 

3 


g 

g 

3 

K 



0 1 g 1 e 1 g 

1 i 1 


-g i 0 i G i g 
( 1 1 


o 

1 1 

<u 
1 1 


” " 1 “ 1 " 1 * 
-g l-g 1-J 1 0 


(53) 


( 54 ) 
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All of the r sub-blocks of II are symmetric. Each of the off-diagonal 
sub-blocks of n is antisymmetric and the negative of the corresponding sub-block 
in the transpose position of II. Hence, the diagonal sub-blocks of H are zero. 


In a similar fashion we write W in the form 

w = ft + w 


(55) 


Here, ft is a diagonal matrix and ft is a symmetric but singly striped 
off-diagonal matrix. These matrices are given in block diagonal form by 


W 


w 


''-(a 


(56) 


^k 


where w is the kth diagonal submatrix of dimension (^)N by (^)N. The matrix W 
has the striped form 


« -2 




0 w 

Z,2Z N 




-w 0 









N 


N 


N 

X 




N 





N 

0 



-w 


(57) 


^k T- 2 

where w is antisymmetric. The form [D WD] leads to r mN by mN submatrix 

T-- 

blocks that fill the entire mrN by mrN array. By Eq. (32), [D WD] has the form 




D^(w' + .,.+w'"’)D 


- D^w’ D 


T~2 
D W D 


D 


■ d^w^-’d 


d^w^-'d 




(58) 


Similarly, [® W2>] has the form 


12) w 2)] = 


0 ^ d'^m^d 
-d'^S^d " 0 


-dVd 


dV-^d 


^0^ D^8(pD 
-o'^8(pD'' 0 


(59) 
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where the matrix M is a combination of the w 's and hence antisymmetric. 

The matrix II as given in Eq. (47) is easily obtained from the numerical matrix 
operations indicated by Eq. (13) with N=l. Upon reformatting IT according to Eqs. 
(50) and (54), and equating these expressions to Eqs. (58) and (59), one 
obtains the appropriate values for the elements of W. For example, equating the 
upper right hand corner sub-matrix block in the two expressions for II , one 
obtains 




(60) 


In this expression § is the synimetric mxm matrix occupying the upper right 
corner block of II. From the definitions of D and W given by Eq. (31) and Eq. 
(56), Eq. (60) may be written as 


eh 


m-l' 




m-l 
.-jr-1 


-W^”^(W^”^+^+» )***-W^' 

1 ^1 m 2m- 3^ 2m- 3 



m-l 


_jjr-l 

2m- 3 



(61) 


‘N^r "• 1 

From Eq. (61) one obtains numerical values for the elements of W for 
arbitrary forms for A^. Analogous equating operations yields the antisymmetric 
elements of W. In practice, the standard deviations of the measurement errors 
usually are uniform over the data streams. Nevertheless, the method described 
here for generating the weights is straightforward and independent of N; the 
computation of II needs to be made only once when m,r,and are specified. 


In summary, W can not in general be a diagonal array, although, it can be 
written as the sum of a diagonal array and a singly striped symmetric array with 
antisymmetric submatrices. The question arises whether or not the least squares 
weighting matrix for the minimally redundant operator given in Eq. (29) also has 
an equivalent diagonal form. However, for this to be true the rows of Eq. (29) 
would have to be the eigenvectors associated with the non-zero eigenvalues of 
Ag» which is generally not true. 
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VII. THE CASE OF OUTAGES 


In the usual tracking situation, satellites rise and set and the mix of 
satellites and terminals in the data streams changes with time. Also, receivers 
occasionally lose data because of instrumental or environmental reasons. These 
cause a considerable sorting problem for the double differencing approach. One 
approach is to decompose the tracking session into subsessions with epoch 
boundaries marked by the onset or termination of an outage for a particular 
terminal. These subsessions would be mutually exclusive and exhaustive. 
Enumerating these subsessions by the index i, it is clear (provided that e is at 
least white across all subsession boundaries) that the covariance matrix for the 
entire tracking session would be given by 

-1 JL -1 

X • X 

1-1 

where A^(i) is the covariance for the ith subsession and n is the total number 
of subsessions. By this means, whatever outages there are among the mr data 
streams will remain invariant during each subsession. 


The outage case for the jth terminal observing the v th satellite can be treated 
by setting o^ = <» in A^. For the case of unifoirm meaurement error statistics, 
other than the ones corresponding to the missing data streams, analytic diagonal 
forms for W can be found, although these expressions become rather complicated 
as the topology of A^ becomes more complex. When the number of outages 
increases sufficiently to reduce a particular terminal down to observations of 
only one satellite (or observations of a particular satellite by only one 
terminal), then the form of W eliminates that terminal (satellite) from the 
solution set and the problem is effectively reduced in dimension from m to m-1 
(or r to r-1). This is, of course, expected from our initial assumptions about 
the stochastic properties of the clocks . 


Some examples of W for outage cases are given below. For the case of a single 
satellite not observed by n on different terminals we set 
and assume a uniform standard deviation of a for the remainder of the data 
streams. A diagonal weighting matrix is obtained (see Appendix 4) in the form 
given in Eg. (56) by the following expressions: 
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w 


V 


1 , 


, r-1 


= - ( 63 ) 

(W + I/m)/(r-l), V = (|) 


where W, a ( 2 >N by ( 2 >N diagonal matrix, is given by 



Here, a,b, and c are given by 


(rnX.^+(r-l+\^)(m-n) ](r-l+\^) 

b = ^ (65) 

rn\^ +(r-l+\^)(m-n) 

c = r-1 + 

nr^X.^+(r-l+\^) (m-n)r 
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Setting X = 0 in these expressions yields the weights for the case of n outages 
involving a single satellite and no outages on the r-1 remaining satellites. 

In the case of two satellites with outages one must consider whether these 
outages involve common terminals, different terminals, or some combination. The 
forms for W are substantially different for the three cases. The diagonal 
weights for these cases are provided in Appendix 4. In general, as the number 
of terminals experiencing outages increases, the topology of W and II becomes 
increasingly more complex resulting in very complicated algebraic expressions 
for W. 

VIII. CONCLUSIONS 


The reader has no doubt detected a preference for linear combination methods 
such as the Householder method rather than double differencing. In the planned 
global GPS tracking network for a demonstration on NASA's planned Ocean 
Topography Experiment mission, TOPEX, [7], 18 GPS satellites would be aloft and 
roughly a dozen globally distributed GPS terminals, including one on TOPEX, 
would be tracking, although not all simultaneously. Nevertheless, ^2)(2) ^ ® 

potentially explosive number. On the other hand, double differencing has a 
number of attractive features that have been noted. It is of interest to 
compare the relative computational times of the two approaches. 

Using the Householder approach it may be shown [8] that the computational time 
T , required to transform Eq. (1) into Eq. (7) at a single epoch (N=l) is given 

n 

approximately by the expression 

Tj^ = H(mr)(m+r+p+l)^ (66) 

where n is a constant coefficient for a specified computer of fixed 
configuration; it accounts for the detailed arithmetic and normalization 
operations in a floating point operation and the transfer of data in and out of 
memory. In the case of fully redundant double differencing one must first form 
the weighted double differenced regression equations from the undifferenced 
regression system and then apply a filter routine. For ease of comparison of 
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the computation times we assume that the Householder transformations are used on 
the weighted double difference array to obtain the triangular array of Eg. (7). 

In this case the computational time, including the double differencing 
operations, is given by 

" ^(^)(|)(P+l)(^/^+P+l) (67) 

Here, \ is a similar constant and \/pi is of the order of unity. Normally, p 
will be some linear function of m and r. One can compare T and T for various 
combinations of these parameters using Eg. (66) and (67) or using actual 
measurements of computation time. In most cases except for modest values 
(m or r ~4) of these parameters the straight Householder approach holds a 
definite computational advantage, typically amounting to a factor of three to 
five. However, the actual computational time for these operations using either 
method is not usually a major concern and is a small fraction of the 
computational time used in data management, editing and validation activities, 
and in forming the regression system itself. Of more concern is the size of the 
arrays that occur in the double differencing approach when the fully redundant 
differencing approach is used. 

An alternative approach is to use the algebraic form for W (when available) to 

T 

generate the matrix operator n. In this case IT is given by |j3> W2>] instead of by 
[3>h ^ ]^3>] thereby avoiding the numerical inversion of the double differenced 
measurement error covariance that results from using a less than fully redundant 

double difference operator. The computation time for the former approach 

2 3 

increases as (mr) , whereas it increases at a rate somewhat higher than (mr) 

for the latter approach. However, when II is explicitly used to operate on the 

undifferenced regression system, one dispenses with explicit double differencing 

in favor of an orthogonal linear combination approach but one saves 

computational time and avoids the large arrays. 

Linear combination methods effected through the Householder or similar 
transformations seem to offer a more attractive approach because they deal with 
only mrN regression rows and readily include the outage cases, and because they 
allow for ancillary stochastic modeling of the clock processes. Both the 
Householder approach and double differencing approaches are currently being 
followed at JPL [ 5 ] . 
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If fully redundant explicit double differencing is the chosen approach, the 
following summarizes the procedures (using a batch mode) to obtain a state 
vector estimate and its covariance, which are equivalent to those obtained from 
linear combination methods: 

(1) Reformat the system of regression equations according to Eg. (3). 

(2) Divide the tracking session into a mutually exclusive and exhaustive 
set of subsessions with the outages remaining invariant within each 
subsession. Each subsession would be characterized by an index i and 
by: m^, the number of participating terminals; r^, the number of 
participating satellites; and N^, the number of simultaneous 
observations by each participant. 

(3) For each of the n subsessions generate the appropriate weighting 
matrix W depending on the form of A^, m^ and r^. 

(4) For each subsession operate on the corresponding regression system 

with ® to generate an expanded but double differenced set of the 
arrays regression equations and form and W^0>y]. 

(5) ’"n the batch mode assemble the weighted and double differenced 
regression systems for all n subsessions into a single ensemble. 

(6) Invert this ensemble to obtain the least squares estimate of x and its 
covariance using one's preferred inversion method such as Householder, 
Cholesky, least squares, etc. 
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APPENDIX 1: ANALYTIC FORMS FOR A"~ FOR THE MINIMALLY REDUNDANT DOUBLE 

DIFFERENCE OPERATOR 

Analytic forms for are available for the minimally redundant double 
difference operator when is factorable into the form given by Eq. (43) or 

(44). In the case where is the minimally redundant double difference 

m — 1 

operator defined in Eqs. (29) and (30), it may be shown that [S>A^ 1 is given 
by 


-1 


T 

[S*Ag2> ] = 


-1 


Y 

0 

• 


r 

r 


■ Y • • • Y ‘ 

• 

• 

• - 1 

0 Y/b^ 

J 

- 

_V = 1 


... 

. r. r J 


(A-1) 


where the matrix Y is defined by 


Y = [DaD ] 


(A-2) 


and where Y ^ is given by 








1 • • • 

1 


0 




;a a 



• 



m 


11 

• 

a - d 

1 m 

Y ^ 

• 








0 

• 

1/a 


4-^ a. 

Ll=l 3j 


# 

1 • • • 

1 



m 




L a a, 
m 1 

a a 
m m 

uniform case (a. = 
3 

= 1), 

Eqn . ( A- 1 ) 

reduces to 


rn 


(r-l)X 

- X 

••• "X 





1 ’■= — 
mr 

-X 

• 

• 

• 

• 






« 


• • 






— X • • 

• 

(r-l)X J 




(A-3) 


(A-4) 


where X is given by 


X = ml -1 ,1 , 

m-1 m-1 m-1 


(A-5) 


and where I and 1 are defined by Eqs. (5) and (6) 
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APPENDIX 2: PROPERTIES OF THE FULLY REDUNDANT DOUBLE DIFFERENCE OPERATOR 


The single difference matrix operator, D^, defined in Eq. (31), can be 
alternatively defined by the recursive relationship 


m+1 


1 I -I 

m I m 

0 1 D 


(A-6) 


with D, = 0 and where 1 and I are defined by Eqs. (5) and (6). The rank of D 
1 mm I '1 \ / V / 

is (m-l)N. It follows from Eq. (A-6) that 

.T 


D^D = ml - 1 r 
mm m mm 


(A-7) 


which can be shown to have( m-l)N identical eigenvalues of value m and N null 
eigenvalues. Using the property 


D 1 z 0 
m m 


(A-8) 


it follows from Eq. (A-7) that the pseudoinverse of D^ is given by 


d'*’ = i D*^ 

m mm 


Here, the pseudoinverse of M has the standard definition 


(A-9) 


M^M 


M M 

= M 

M 

t T T 

= (M^) M 

M^M 

= M^(M^)^ 


(A-IO) 
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A property of D for n < m, given by the relationship 


0 , 0 


mn 0 (dTd 

* m m 


0 , 0 


0 pTo 

’ n n 


0 I 0 


O P^D 
* n n 


is useful in "spectral decomposition" analyses; for example, a sub-matrix block 
of n in certain cases may be expressed as 


A = > a 




0 r m-k 

dT 1 k 


The double difference operator 2>is defined by Eq. (32). The inner product 
T 

Q), is given by 


2)'0i = r 


D D 0 

\ 

N 

\ 

\ 

s 

S 

0 D^D 


D^D- 


- D^D 


(A-13; 


or alternatively, using Eg. (A-7), it may be expressed as 


Qi Q) = m r I -r 
mr 


11' 0 
m m 


0 1 r 

m m 


- m 

I 0 

m 

V 

V 

\ 

+ 


1 

o 

1 



mm mm 

1 ’i^ 1 i' 

mm mm 


It follows that the rank of® ®is (m-l)(r-l)N. From Eqs. (A-9), (31), and (A-13) 
it follows that the pseudoinverse of ® is given by 


t T 

S) = 

mr 


(A-15) 


Also,2)'^/mr has (m-l)(r-l)N unit eigenvalues and (m+r-l)N null eigenvalues. 
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APPENDIX 3: THE STRUCTURE OF THE H -MATRIX FOR A UNIFORM WEIGHTING MATRIX 

To obtain insight into the structure of II in double differencing, it is helpful 

to examine the structure of n that follows from its direct computation using Eg. 

(13) for the case of uniform measurement error covariance. To examine II 

further, we will need the pseudoinverse of [B B]. Alternately, one could 

strike one of the columns of B, thereby eliminating the rank deficiency in B, 

T -1 

and deal with the inverse of [B B] instead of the pseudo inverse . Column 

striking spoils the symmetry of the problem as posed, although it leads to the 

same results as the pseudoinverse applied to the fully dimensioned B. Either 
T -1 

way, [B A B] is difficult to invert except in special cases where A is 

® T -1 ® 

factorable, outages, etc. Let us show the form of [B A B] when A is 

6 © 

uniform. It can be shown for this case that 


• o* 


r 1 

!-i 

m 

1 m r 

— 

1 

-<f 

-1 

! ml 

f m 

1 r 


1 


(A-16) 


and that its pseudoinverse is given by 




mr (m+f) 








j r(m+r)*I^-rfr+am) 


(A-17) 


After some matrix algebra it may be shown using Eq. (13) that II becomes 










j — j 

• • • 

1 

1 




— 1 1— 

1 1 

• • • 

1 

H- 

- ml + 1 1^ 

m fii 

• 


1 * 1 

• 

1 

1 

• 

• 


1 • j 
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1 

1 

• 
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1 • 1 

• 

1 

1 

• 



1 -"*> ♦' 1 
1 Mn ifl Ifl 1 

• • • 







J 


(A- 18) 


But, using Eq. (A-7) we now identify the array in Eq. (A-18) with that in Eq. 
(A- 13). In terms of the 9>-matrices , we conclude that II is given by 
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(A-19) 


2 _ 

n - ^ Q^Q 
mr 

Thus, the double differencing structure appears in the expanded state space 
array for the covariance of x when it is factored according to Eg. (A-19). 

One can also derive other analytical forms for W from the expanded state space 
approach whenever explicit forms for the submatrix elements of IT can be obtained 
analytically from Eg. (13). Another approach using the idempotent properties of 
IIa| , is discussed in Appendix 4. 


30 



APPENDIX 4: OTHER DIAGONAL WEIGHTING MATRICES FOR MULTIPLE OUTAGES 


Inasmuch as the major topological properties of H are invariant to N, we set N=1 
in the discussion to follow without significant loss of generality. Another 
useful way to arrive at explicit forms for II is to use the relationship 


n = n Ag n 

which follows from Eq. (13); also 

Trace n ] = (m-l)(r-l) 


(A-20) 


(A-21) 


2 2 2 2 

In addition, H satisfies the m r -(ra-1) (r-1) 
by 

ne = b'^ii s 0 


orthogonality conditions defined 

(A-22) 


In theory these conditions are sufficient to solve for II when the form for 

A is given. In practice, the topology of II, except for the simplest of forms 

for Ag, becomes so complex that the resulting explicit algebraic expressions 

for the elements of H in terms of m,r and the elements of A are rather tedious 

e 

to evaluate. In these cases a more practical approach is to numerically invert 
n for N=1 as given by Eq. (13). 

As examples of this analytical technique the single and double outage cases are 
presented. For a single satellite outage we consider a A^ with the form 



(A-23) 
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where X = 0 corresponds to the case of n terminals not observing a single 
particular satellite with the remainder of the terminals observing all r 
satellites. Applying part of the conditions in Eg. (A-22) one obtains a IT of the 
form 


n = 


(r-l)A -A - A 

-A A+(r-2)B -B . . . B 

• -B • A+(r-2)B 


L -A -B A+(r-2)B •* 

Applying Eg. (A-20) one obtains the conditions on A and B 


(A-24) 


rA^ = A 


(A-25) 


and 


(r-l)B = [D D/m - A] 


(A-26) 


Next, eguating the elements of the matrices in these expressions and applying 
Eg. (A-21) and the remainder of the orthogonality conditions one obtains for A 
the form 


A = 



m 



b 

to 

T 

b 11 

T 

dl+cll 




- 


n 

m-n 


(A-27) 


where the coefficients a,b, and c are given in Eg. (65), and d and r are given 
by. 


e = X^/(r-l+X^) 
d = 1/r 


(A-28) 


The elements of B follow from Eq. (A-26). Note that no syinmetricizing of A and 
B is necessary because they are already symmetric for this special case. With 
these expressions for A and B placed in Jl as given by Eq. (A-24) which is in 
turn eguated to [S> W^] as given by Egs. (58) and (61), one obtains the algebraic 
expressions for the diagonal weights that are given in Eg. (65). 

For the double outage case we let one satellite incur outages on x different 
terminals. We let a second satellite incur y outages that involve terminals 
among those x terminals not observing the first satellite. We let the second 
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satellite also incur z outages involving terminals observing the first 
satellite. Thus, the second satellite incurs a total of y+z outages with y<=x 
and z<=m-x. (There is a duality symmetry in this problem that allows one to 
exchange the roles of satellites and terminals.) The form for may be given 
by 


-1 

Ae 



m-x+m(r-2) 

, y 

m-y-z 


(A-29) 


Following the same routine as used in the previous example, it may be shown that 
the topology of II has the form 
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1 
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“Cl 

- E^ 


\ 
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“E| 

A+(r-2)E^ 

• 


1 

1 

• 


(A-30) 


Where each of the sub-matrices A,B,C and E are m by m. Here, the x missing 
data streams for the first satellite correspond to zeros in the first x rows and 
columns of II . Similarly, the z missing data streams for the second satellite 
(which for convenience has been placed last) correspond to zeros in the last z 
rows and columns. The y outages of the second satellite correspond to zeros in 
rows and columns numbering (r-l)m+l through (r-l)m+y. Thus B and E, the top and 
bottom edge matrices are of rank m-x-1 and m-y-z-1, respectively. The rank of A, 
the corner matrix, is the smaller of the ranks of A and E. These matrices are 
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generally not symmetric and will have to be decomposed into symmetric and 
antisymmetric parts before equating them to the diagonal weight matrices given 
in Eqs. (58) and (59). The core matrix C, is symmetric and explicitly appears 
whenever r>=4. For r=2 or 3, C may be considered as a virtual matrix which does 
not explicitly appear in n . For ail values of r^2 it may be shown upon 
applying Eqn. (A-20) that B,C and E satisfy the relation 

(r-2)C = D^D/m - B-E (A-31) 

The topologies of B,C and E can be described in block matrix form by 


0 
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0 

0 
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0 
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0 
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(A-32) 


(A-33) 


(A-34) 
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It may be shown that the non-zero elements of B satisfy the relation 


b. . 
ID 


c . . 
ID 



(A-35) 


Similarly, the non-zero elements of E are given by 

'i3 ' 'ij * [e *E 'kj 

k=l k=m--z+l 

Also, the non-zero elements of A are given by 

m 

a.. = b. . + b., /(m-y-z) 

1] / j ik ' ' 

K=in-z+l 

It may be shown that the form of C has the additional properties given by 


(A-36) 


(A-37) 
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T 

bll 

T 
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cl+dll" 

T 
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T 

bll 

T 
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fll 

T 

jll 

ki+m" 
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Thus, it suffices to determine the elements of the core matrix C given in Eq. 
(A-38). Applying these relationships between B,C and E in Eqs. (A-31), (A-35) 
and (A-36) and the remainder of the orthogonality conditions in Eq. (A-22), one 
obtains 


a = l/(r-2) 

c = l/(r-l) 

g = 1/r 

k = l/(r-l) 

b = -l/m(r-2) 

d = -1 (r-1) (v+z) (m-x) - z(m-x+y) 

m(r-l) m (r-l)(r-2)H 


V (A-39) 


e = (r-2 ) (m-x) (rz-(r-l )m) + y[(r-l) (m-x)-zl 

m (r-2)H 
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f 


-r(r-2)(m-y-z)(m-x) + y(m-x-2) 
m(r-2)H 


h = - (r-2)[(r-l)(m-x)-x][(r-l)(m-z)- )1 + y[ (r-1) (r(r-l) (m-x)-(r-2)m)-rz] 

mr(r-2)H 


j = (r-2) (iti-z) f (r-1) (m-x)-xl + yf(r-l) (m-x)-(r-2)m-z1 

m(r-2)H 


i ~ - 1 + rx(m- 2 ) 

m(r-l) m(r-l)H 


+ yfm-z-(r-l) x] 
m(r-l) (r-2 )H 


where 


H 


(r-l)‘ 


(m-x) (m-y-z)-z(x-Y) 


(A-40) 


The expressions for the elements of A,B and E follow from Eqs. (A-35), (A-36) 
and (A-37). 


As an example, the matrix II, obtained numerically from Eq. (13), is presented 
below for the case where m=5, r=4,x=2,y=l, and z = 1. 


80n for m,r,x,y,z = 5,4,2, 1,1 
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0 0 0 0 0 

0 

0 0 0 

0 

0 

0 

0 0 

0 

0 0 

0 

0 0 

0 0 0 0 0 

0 

0 0 0 

0 

0 

0 

39-21- 

18 

0 -1- 

14 

6 9 

0 -1-14 6 9 

0 

2-11 9 

0 

0 

0- 

21 39- 

18 

0 -1 

6- 

14 9 

0 -1 6-14 9 

0 

2 9-11 

0 


0- 

18-18 

36 

0 2 

8 

8-18 

028 8-18 

0 

-4 2 2 

0 

0 

0 

0 0 

0 

32 -8 

-8 

-8 -8- 

-32 8 8 8 8 

0 

0 0 0 

0 

0 

0 

-1 -1 

2 

-8 41- 

12- 

12 -9 

8-23 447 

0- 

18 9 9 

0 

0 

0- 

14 6 

8 

-8-12 

46- 

14-12 

8 4-18 2 4 

0 

8-14 6 

0 

0 

0 

6-14 

8 

-8-12- 

14 

46-12 

8 4 2-18 4 

0 

8 6-14 

0 

0 

0 

9 9- 

18 

-8 -9- 

12- 

12 41 

874 4-23 

0 

2 -1 -1 

0 

0 

0 

0 0 

0- 

-32 8 

8 

8 8 

32 -8 -8 -8 -8 

0 

0 0 0 

0 

0 

0 

-1 -1 

2 

8-23 

4 

4 7 

-8 41-12-12 -9 

0- 

18 9 9 

0 

0 

0- 

14 6 

8 

8 4- 

18 

2 4 

-8-12 46-14-12 

0 

8-14 6 

0 

0 

0 

6-14 

8 

8 4 

2- 

18 4 

-8-12-14 46-12 

0 

8 6-14 

0 

0 

0 

9 9- 

18 

8 7 

4 

4-23 

-8 -9-12-12 41 

0 

2 -1 -1 

0 

0 

0 

0 0 

0 

0 0 

0 

0 0 

0 0 0 0 0 

0 

0 0 0 

0 

0 

0 

2 2 

-4 

0-18 

8 

8 2 

0-18 882 

0 

36-18-18 

0 

0 

0- 

11 9 

2 

0 9- 

14 

6 -1 

0 9-14 6 -1 

0- 

18 39-21 

0 

0 

0 

9-11 

2 

0 9 

6- 

14 -1 

0 9 6-14 -1 

0- 

18-21 39 

0 

0 

0 

0 0 

0 

0 0 

0 

0 0 

0 0 0 0 0 

' 0 

0 0 0 

0 
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